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Approximate  equations  for  a Boussinesq  model  with  viscous  dissipation  and 
thermal  conduction  describing  buoyant  convection  driven  by  a heat  source  in  a 
rectangular  enclosure  are  derived.  The  finite-difference  algorithm  for  computing 
transient  solutions  in  two  dimensions  to  these  equations  is  presented.  The  al- 
gorithm allows  the  enclosure  fluid  to  be  stratified  in  a direction  parallel  to  the 
enclosure  walls  initially,  or  for  gravity  to  have  an  arbitrary  direction  relative  to 
the  enclosure  (but  with  no  initial  stratification).  Computational  results  of  tran- 
sient, two-dimensional  buoyant  convection  for  very  high  resolution  are  presented. 
The  hydrodynamics  is  directly  based  on  the  time-dependent  Navier-Stokes  equa- 
tions; the  model  is  valid  in  the  Boussinesq  approximation.  No  turbulence  model 
or  other  empirical  parameters  are  introduced.  There  is  no  inflow  or  outflow  at 
boundaries;  this  assumption,  although  rather  restrictive,  allows  the  mathematical 
problem  to  be  properly  formulated  so  that  no  other  empiricism  is  introduced  by 
specification  of  the  algorithmic  boundary  conditions.  A finite-difference  scheme 
second-order  in  space  and  first-order  in  time  is  used  to  integrate  the  evolution 
equations,  and  an  elliptic  solver  is  used  to  solve  the  pressure  equation.  The  algo- 
rithms have  been  verified  by  comparisons  with  exact  solutions  to  the  equations  in 
simple,  special  cases,  and  predictions  of  the  overall  model  when  the  viscosity  and 
thermal  conductivity  are  zero  have  been  compared  with  experimental  results.  The 
use  of  Lagrangian  particle  tracking  allows  one  to  visualize  the  flow  patterns. 


1 Introduction 

The  authors  have  pubhshed  previously  a description  of  the  mathematical 
model  and  their  algorithm  for  computation  of  the  buoyant  convection  in- 
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duced  by  a fire  evolving  in  a room  [1]  - [7].  The  model  is  for  a dissipation 
free,  thermally  expandable  fluid,  i.e.,  one  for  which  density  and  tempera- 
ture variations  can  be  large,  but  pressure  variations  are  small  [1],  and  for 
a Boussinesq  model,  one  in  which  the  density  variations  are  also  small. 
The  algorithm  can  be  generalized  easily  to  include  viscous  dissipation  and 
thermal  conductivity.  Previously  the  authors  had  argued  that  one  could 
not  resolve  both  large-scale  motions  associated  with  room-scale  buoyant 
convection  and  motions  at  the  dissipative  length  scale  with  computational 
resources  available.  However,  with  the  revolution  in  computational  capa- 
bihty  currently  taking  place,  this  statement  is  no  longer  strictly  true,  at 
lesLst  in  two  dimensions. 

In  this  document  the  authors  generalize  their  model  and  algorithm  in 
two  dimensions  to  include  viscous  dissipation  and  thermal  conduction.  Also 
included  are  recently  introduced  generalizations  which  allow  for  either  an 
initial  background  stratification  in  a direction  parallel  to  the  enclosure  walls 
or  for  the  enclosure  to  have  an  arbitrary  direction  relative  to  gravity  (but 
not  both  in  the  same  calculation).  We  present  the  details  of  the  algo- 
rithm in  this  report  (which  may  be  regarded  a^  a working  document  on 
the  two-dimensional  boussinesq  code  with  dissipation).  Because  the  model 
is  restricted  to  two  dimensions,  very  high  resolution  computations  can  be 
performed,  which  allow  us  to  resolve  both  large-scale  buoyant  convection 
and  small  scale  dissipation  for  Reynolds  numbers  of  interest  for  enclosure 
fires.  The  model  of  buoyant  convection  is  presented  in  Section  2 and  a brief 
description  of  the  numerical  methods  in  Section  3.  Section  4 contains  a 
detailed  description  of  the  algorithm.  Some  sample  results,  together  with 
discussion  of  the  physical  phenomena  observed,  are  described  in  Section 
5.  An  Appendix  is  also  included  which  contains  the  algorithm  used  to 
introduce  and  track  the  Lagrangian  particles. 


2 Hydrodynamic  Model 

Traditionally,  two  approaches  to  the  computation  of  fire-induced  buoy- 
ant convection  have  been  reported:  direct  integration  of  the  Navier-Stokes 
equations  using  molecular  values  for  viscosity  and  thermal  conductivity  or 
integration  of  these  equations  using  a turbulent  viscosity  and  conductiv- 
ity to  account  for  fluctuations  occurring  at  the  large  Reynolds  numbers  of 
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practical  interest.  The  former  approach,  although  most  desirable,  is  not 
practically  feasible  in  three  dimensions  with  today’s  computers;  computers 
are  not  big  and  fast  enough  to  calculate  simultaneously  the  largest  scales 
of  buoyant  convection  and  the  smallest  scales  at  which  dissipation  occurs. 
Alternately,  the  use  of  a turbulence  model  in  the  equations  introduces  func- 
tional forms  and  empirical  constants  which  do  not  have  a fundamental  the- 
oretical basis  at  this  time.  However,  in  two  dimensions,  direct  simulation 
of  the  Navier-Stokes  equations  at  Reynolds  numbers  of  practical  interest 
for  fire-driven  flows  is  possible  and  even  practical,  and  these  simulations 
are  the  subject  of  this  report. 

The  essence  of  the  buoyant  convection  model  is  described  as  follows. 
We  consider  a Boussinesq  fluid  with  constant  \’iscosity  and  thermal  con- 
ductivity in  a rectangidar  enclosure  driven  by  a prescribed  heat  source.  We 
start  with  the  equations  of  motion  for  a thermally  expandable  ideal  gas 
[1]  in  which  we  include  viscous  dissipation  and  thermal  conductivity  (with 
constant  viscous  and  conduction  coefficients). 


5 , . 

/ du,  dui\  dp 
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d^Ui 
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dT 
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= Q + K 


d^T 
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Po  = pRT 


(1) 


Here,  all  symbols  have  their  usual  fluid  dynamical  meaning:  p is  density,  u , 
are  the  velocity  components,  p is  pressure,  g is  the  acceleration  of  gravity, 
ki  are  the  components  of  the  vector  describing  the  direction  of  the  gravity 
vector,  u is  kinematic  viscosity,  Cp  is  the  constant-pressure  specific  heat,  T 
is  temperature  and  K is  the  thermal  conductivity,  t is  time  and  Q is  the 
spatially  and  temporarly  prescribed  heat  source. 

If  we  combine  these  equations  as  described  in  [2],  nondimensionahze 
them  a.s  described  in  the  appendix  of  this  reference  and  make  the  Boussinesq 
approximation,  we  obtain  the  following  set  of  equations. 
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dui  dui 

w + + 


dp 

dxi 


dp  dp  _ 

di  dxi 

{p-po{z))ki  = 

dui  _ 
dxi 


-^!—^Q  + kV^P 
7 

vV^Ui 

0 


(2) 


Here,  all  symbols  have  the  meanings  given  above  but  in  dimensionless  form, 
Pq{z)  is  the  initially  stratified  density  profile  assumed  to  depend  only  upon 
the  vertical  coordinate  2r,  and  k is  the  dimensionless  thermal  conductivity. 
The  dimensionless  quantities  are  defined  as  follows:  lengths  are  relative  to 
the  height  of  the  enclosure,  pressure  is  relative  to  ambient  pressure,  time  is 
relative  to  the  height  of  the  enclosure  and  the  velocity  scale,  and  velocity 
is  relative  to  a scale  U,  and  the  density  perturbation  is  relative  to  ambient 
density  with  a small  parameter  (3,  If  we  denote  temporarily  the  dimensional 
quantities  by  an  asterisk,  then  these  relationships  can  be  written  as  follows: 


X*  = Hx„  r = 

u*  = U Ui 
p-  = po(0)J7V 
p-  = po(0)(l  + I3p) 

Po(^')  = po(0)(l  + 0p(z)) 


U and  /3  are  defined  in  terms  of  the  magnitude  of  the  heat  source  as  follows: 

/3  = uy{gH) 

U = {Qogl{poC,ToH)Y/^  (3) 

U = {qog/ipoC,To)y>^ 

Here  Qq  is  the  strength  of  the  three-dimensional  heat  source  in  units  of 
energy  per  unit  time,  qo  is  the  strength  of  the  two-dimensional  heat  source 
in  energy  per  unit  length  per  unit  time,  H is  the  height  of  the  enclosure 
and  all  other  quantities  have  their  usual  meanings.  See  [2]  and  [6]  for 
more  details  of  the  scaling.  Henceforth,  all  variables  will  be  regarded  as 
dimensionless. 
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In  the  two-dimensional  case,  these  equations  can  be  rewritten  as  follows: 


dp 

Jr  F F 

^ PX  ^ PZ 


du 

di 


dp 


+ + T P9x 


dw  ^ dp  ^ 
du  dw 


/ d^u  \ 
dz'^ ) 

_ / d^w  d^w\ 


where 


P = P-  Po(2) 

ax 

/ dpo  dp\ 


Fi  = (1/2)^  - wu 

F,  = (l/2)^  + «u, 

dw  du 
dx  dz 
= u^  + w^ 


(4) 


Boundary  conditions  for  these  equations  are  that  there  be  no  inflow 
or  outflow  at  boimdaries,  that  there  be  a no-slip  or  free-slip  conditions 
at  boundary  walls,  and  that  the  walls  are  adiabatic  or  kept  at  a constant 
temperature.  The  initial  conditions  are  that  the  fluid  is  quiescent. 


3 Numerical  Methods 

Equations  (1)  are  a mixed  parabolic/elliptic  system  of  partial  differential 
equations;  i.e.,  the  equations  for  the  density  and  for  the  velocity  components 
are  parabolic,  whereas  that  for  the  pressure  is  elliptic.  The  incompressible 
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equations  of  hydrodynamics  are  well  known  to  have  this  mixed  character. 
When  there  is  no  dissipation  in  these  equations,  it  is  important  not  to  intro- 
duce any  through  the  numerical  scheme.  Analytical  studies  of  the  abihty 
of  several  candidate  finite  difference  schemes  to  calculate  internal  gravity 
waves  [3]  led  to  the  conclusion  that  methods  of  second  order  accuracy  in 
space  and  time  would  be  necessary;  without  dissipation,  the  scheme  cho- 
sen is  dispersive,  but  not  dissipative.  All  time  derivatives  are  replaced  by 
central  differences  over  twice  the  time  step  size  (a  leap-frog  scheme)  while 
all  viscous  and  conduction  terms  are  either  differenced  at  the  lagged  time 
level  in  one  case,  or  according  to  a DuFort-Frankel  scheme  in  the  other. 
Both  schemes,  the  lagged-dissipation  one  and  the  DuFort-Frankel,  have 
been  used.  Other  terms  in  the  evolution  equations  (the  first  two  of  Eqs. 
(1))  are  in  general  evaluated  at  the  mid  level  of  the  three  level  scheme. 

The  spatial  grid  is  taken  to  be  uniform  in  each  of  the  two  directions, 
although  the  mesh  length  may  be  different  in  each  direction.  Within  each 
rectangular  mesh  cell,  vector  components  are  evaluated  at  the  sides  and 
scalar  quantities  at  the  center  of  the  cell.  The  staggered  grid  permits  central 
differences  to  second  order  accuracy  for  all  lineax  operations.  The  nonlinear 
terms  must  be  considered  separately.  The  density  evolution  equation  in 
continuous  form  is  the  maiss  conservation  equation  minus  the  expression 
for  the  velocity  divergence.  Each  of  these  two  equations  is  approximated 
by  central  differrences  and  then  subtracted.  The  density  at  all  faces  is 
approximated  by  the  mean  of  the  density  at  the  centers  of  adjacent  cells. 
This  procedure  ensures  global  mass  conservation  ais  well  as  second  order 
accuracy. 

The  momentum  equation  is  differenced  in  the  vector  invariant  form 
shown  in  Eqs.  (1).  This  ensures  nonlinear  stability  and  complete  compat- 
ibility between  the  “primitive  variable”  formulation  presented  here  and  a 
vorticity,  stream-function  formulation  (in  the  dissipation- free  c2Lse),  see  [4] 
for  details. 

The  pressure  equation  is  the  discretized  version  of  the  time  derivative  of 
the  third  of  Eqs.  (1),  using  the  central  difference  approximation  to  the  di- 
vergence of  the  velocity  and  with  the  time  difference  of  the  velocity  replaced 
using  the  discretized  momentum  equations.  Mathematically,  the  calcula- 
tion of  the  pressure  requires  the  solution  of  an  elliptic  partial  differential 
equation.  For  the  Boussinesq  model,  the  Hnear  algebraic  system  arising 
from  its  discretization  has  constant  coefficients  and  can  be  solved  by  a fast 


6 


direct  method,  see  [6]  for  details.  The  solution  to  the  pressure  equation 
constitutes  the  bulk  of  the  numerical  computation  since  the  density  and 
the  velocity  are  updated  expHcitly  once  the  pressure  gradients  are  known. 

Finally,  stabihty  of  the  computational  scheme  imposes  a limit  on  the 
time  step  size  relative  to  the  spatial  mesh  sizes,  see  [2]  and  [3].  Also,  ac- 
curacy of  the  computation  is  an  important  consideration,  which  has  been 
examined  in  [2],  [3]  and  [6],  and  verification  that  the  numerical  methods 
solve  the  partial  differential  equations,  at  least  in  special  cases,  has  been 
addressed  in  these  references.  In  addition,  the  basic  features,  such  as  the 
plume  rise- time  in  a uniform  density  environment,  have  been  compared 
with  experimented  results  to  verify  the  basic  predictive  capability  of  the 
dissipation- free  model  [5].  With  such  careful  consideration  of  the  basic 
model  and  numerical  methodology,  we  have  confidence  in  the  predictions 
made  by  these  computations,  and  feel  justified  in  interpreting  physical  fea- 
tures arising  in  them. 


4 The  Algorithm 


4.1  The  Density  Equation 


We  now  write  out  the  details  of  the  algorithm.  First,  for  the  density  equa- 
tion, return  to  the  full  continuity  equation  (in  dimensional  form): 


dp  ^ / X ^ 

_+_(p„.)  = o 

dp  dp  dui  _ 


(5) 


Use  the  nondimensionalization  presented  earlier 


dp  djpo  -h  p) 

di  dxi 


+ [1  -h  (3(po  -t- 


= 0 


(6) 


where  p = l-hyd(po“hp),  where  now  all  variables  are  dimensionless  and 
where  we  have  divided  through  hy  /3.  If  we  now  formally  allow  P —*■  0^  then 
we  recover  the  Boussinesq  equation  for  the  density. 


^ d{po  + p)  duj 

di  dxi  dxi 


(7) 
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The  same  procedure  used  above  for  the  partial  differential  equations 
(PDE)  can  now  be  applied  to  the  finite  difference  equations  (FDE).  We 
desire  to  keep  conservation  form  for  the  FDE  as  well  as  for  the  PDE.  Hence, 


dpik  , rpi+l,k  + Pik  Pik  + Pi-l,k 


dt 


+ ( 


‘^ik 


1,A:J  c 

OX 


, rPt,fc+l  + Pik 

+ Wik 


Pik  + Pi,k-\  T 1 p, 

2 “^•■*-157  = 0 


or,  rewriting, 
dpik 


Pi+l,k^ik 


di  2Sx 
, Pi,k+l'Wik 


Pi-l,kUi-l,k  Pik  / X 

+ :rr-[Uxk  - Ui-i,k) 


28x 


28x 


28z 

Using  the  identity 


Pi,k-\'^i,k-\  Pik  / \ n 

+ :ZT-(^ik  - Wi^k-l)  = 0 


28z 


28z 


ac  — hd  = [(a  + h){c  — d)  + (a  — 6)(c  -h  d)]/2 


we  find 


Pi-\,k){'^ik  — y'i-\,k) 

“1“  Pi—l,k^(^^ik  “1“  ^i— l,fc)  2,8x  l,fc) 

4-=^(p,-,*:+l  + Pi,k-l)(Wik  - Wi,k~l) 

^dz 

4- - Pi,k-i){wik  + Wi,k-i)  + “ 'Wi^k-i)  = 0 


or 


where 


dpik 

di 


+ ■F'pxifc  + Fpzik  H"  (^/^)Pik^ik  — 0 


A,jt  = 


^^3  / ik 


7Po  Rp^ 


J ik 


Rpxik  — Pi-1, k)  i^i,k  1,A:) 


(8) 


p "TT  {pi+l,k  Pi—l,k)  {^i,k  4“  1,/:) 

48  X 
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Fpzik  — (P»,fc+1  Pj,A:+i)  ^t.fc— l) 


{PiMl  - Pi,k-l){^i,k  + U?,- Jb_i) 


A8z 


Finally,  use  the  same  nondimensionalization  as  cited  above,  using  p = 
1 4-  d{po,ik  + pik)‘  Note  that  only  pik  and  p.jt  depend  upon  time.  Also,  note 
that,  we  can  divide  by  l3  and  then  allow  it  formally  to  vanish  as  done  above 
for  the  PDE  to  obtain  the  finite  difference  equation  for  the  Boussinesq 
model. 

Fpnk  + Fpzxk  + Dik  = 0 (9) 

where 


Fpxik  — (pO,i+l,Jb  + P0,i-l,jb  4"  Pt+l,/fc  4-  pi-l,k)  {^i,k  ” 

4”  4(5^  (P0,»+l,fc  Po,i—l,k  4"  Pi-\-\,k  Pi—\,k)  {^i,k  4“  ^i—l,k^ 

■F* ozik  ~ ^ (Po,i,A:+l  4"  Po,i,k+l  4“  Pi,k+1  4"  Pi,k—l)  {^i,k  l) 


1 


‘^Sz 


{Po,i,k+\  Po,i,k—l  4“  Pi,k+1  Pi,fc— l)  4”  ^i,k—l^ 


and  where  Dik  is  the  dimensionless  form  of  A,jt  divided  by  /3, 

7-  1 Q 


Dik  = 


7 Po 


- «vV) 


ik 


Here  k = KKpCp). 

At  boundaries,  the  density  fluxes  are  determined  by  the  no-inflow,  no- 
outflow conditions.  Also,  we  must  specify  adiabatic  or  cold-wall  boundary 
conditions,  which  determine  the  temperature  and  hence  the  density  (since 
the  perturbation  density  is  the  negative  of  the  perturbation  temperature) 
in  ficticious  cells  adjacent  to  the  boundaries. 

At  i = 1: 


Dpxlk  = 2^  (P0,2,A:  4-  P2,k)  U2,k 
Po.Jk  = 


(10) 


9 


At  i = I: 


^ oxlk  — 


2Sx 


{pO,I-l,k  + pi-l,k)  Ul-l,k 
pl+l,k  = ±p/,jfc 


(11) 


At  A;  = 1: 


Fpzil  — Pi,2)'^i,2 

Pi,0  = ±Pi,i 


At  k = K: 


^ oziK  — 


^^{pO,i,K-l  + Pi,K-l)^i,K-l 
Pi,K+l  = ^Pi,K 


(12) 


(13) 


Here,  the  plus  sign  corresponds  to  adiabatic  boundary  conditions  (BC)  and 
the  minus  sign  to  cold  wall  BC.  Adiabatic  BC  imply  zero  derivative  in  the 
temperature  (density)  perturbation  across  the  boundary  whereas  cold  wall 
BC  imply  the  temperature  (density)  perturbation  is  zero  at  the  boundary. 

Two  different  discretizations  in  time  have  been  used.  In  each  case,  the 
discretizations  of  the  time  derivative  has  been  chosen  as  a leap  frog. 

^ dt  ’ 2SV^''‘  ' 

In  the  first  case,  we  use  a lagged-diffusion  temporal  discretization  for  the 
dissipation. 


= J^{P0,i,k+l  "*  2/?0,iA:  + P0,i,k-l)  “ + P7k^) 

+ P'i-lk)  + Plk-l) 


Sx^ 


Sz^ 


Define 


^p,tk  — ~^px,ik  ~ ^pz,ik  Q?k  + ’^^{P0,i,k+1  '^p0,ik  + Po,i,k-l) 


px,ik 


pz,ik 


-2(^  + + TZ^ip’i+lk  + P"-u)  + + Plk\) 


8x'^ 


8z^ 
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We  also  define  to  be  the  same  as  but  with  all  quantities  defined 
at  time  level  n.  In  this  case,  the  temporal  discretization  for  second-order 
time  steps  yields 


= 28tRi,,  + 


(14) 


and  for  the  first-order  time  steps  (initially) 


pT  = + Plk 


(15) 


On  the  other  hand,  when  we  use  the  DuFort-Frankel  temporal  dis- 
cretization scheme 


(vv)r*  = (v^w^)+p])”. 

= — 2po,.fc  + PO,t-,fc-l)  - + p'xk^) 


8x^ 

In  this  case,  we  define 


px,ik  pz,ik 


^p,ik  ~ ^px,ik  ^pz,ik  ^ Q ik  2po,tA:  + PO,i,A:-l ) 

(P"+I,jb  + ?i-\,k)  + TZl^^tk^X  + Kk-\) 


Sx^ 


Sz2 


The  temporal  distretization  yields 
1 


We  define 


r,  = 2KSt(±+±) 


Then,  for  the  second-order  time  steps 

P^'  = (1  + r,r^l2StK.,  + (1  - 

and  for  the  first-order  time  steps  (initially) 

p?<t'  = + (1  - ^.)pa 


(16) 


(17) 


11 


In  these  calculations,  the  following  forms  for  the  heat  source  have  been 
used: 

Q?k  = rQr..Q..k  (18) 

where  /"is  either 

/"  = Asinh(ar) 
or 

f'^  = B sin(Q^") 

for  0 < ^ < (ir/a),  and  /"  = 0 for  t > n/ a . Also, 

Qx,x  = \/^7/^exp[-l3j,{xi  - Xc)^] 

Qz,k  = \/l3z/7Texp[-l3^(zk  - Zcf] 

or 

Q^,k  = Aexpf-A^rjt] 

4.2  The  Momentum  Equations 

In  the  velocity  equations,  we  use  the  following  difference  scheme: 
du 

+ Fxik  + {Pi+i,k  - Pi,k)/{Sx)  - gx{pi+i,k  + pi,k)/2  = u{V'^u)ik 
dxv 

+ F,ik  -h  (pi,k+i  - Pi,k)l{^z)  - gz{pi,k+i  + P»-,A:)/2  = v{V'^w)ik  (19) 

The  velocity  fluxes  are  defined  as  follows: 

Fxxk  = + 0.6{w^ikUJik  + w^i,k-i<-^i,k-i) 

F zik  ~ (kQx,k+l  ~ “1“  ^u;*— l,fc) 

where 

qfk  = 0.25[(u,ib  + -h  {wik  + 

u;,A:  = {Ui,k+1  - Uik)l{Sz)  - {Wi^i^k  - w,k)l{Sx) 

Wujik  = 0.5  * {wik  + Wi+i^k) 

Uujik  = 0.5  * {Uik  4-  Ui^k+i) 
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The  temporal  discretization  is  done  exactly  els  for  the  density  equation;  we 
either  used  a lagged  dissipation  or  a DuFort  Frankel  scheme.  In  both  cases, 
the  temporal  derivatives  are  handled  by  a leap-frog  discretization 


1 

28i 


When  we  use  a lagged  dissipation  method,  we  define  the  dissipation  terms 
as  follows: 


6z^ 


Sx^ 


n 

tk 


1 


(V  + ^?-i,k)  + 777(^a+i  + ^i,k-i)  “ 2(  — + 


62^ 


8x^  82^ 


And  we  define  the  right  hand  sides  as  follows: 


u?,  = -f:„  - (pr+:,. 
wr,  = -F^.,  - 


“ P?,k)li^^)  + 9z(p'i+i,k  + pZk)!"^  + 
“ Ptk)/(^^)  + dzip'i.k+l  + p'i.k)!^  + 


urk  = -F"k  - ip?+i,k  - pWKSx)  + g,{p"+,,k  + PlkM'i  + 

WPk  S -F:,k  - {p?Mi  - Plk)/{ix)  + + Plk)l2  + 

+u(V^w)lk 

For  the  second-order  time  steps  using  the  lagged-dissipation  scheme, 

ur+i  = 28tU",  + 

= 28tWl‘^  + U'.T' 
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and  for  the  first-order  time  steps  (initially)  using  the  lagged  dissipation 

= stwr, + w". 


In  the  DuFort-Frankel  scheme, 


+ K-i.k)  + j^^Kk+i  + “.Vi)  + ‘ 


Sz^ 

(V^U))";,  = ^(“’"+1,*  + W^-l,k)  + ^(“'.>+1  + ^?,k-l)  + (^  + 
In  this  case,  we  define 


tk 

n+l 


k 

n-1 


* +«':r‘) 


^ik  — 


wr,= 


~^xik  ~ {P’i+l,k  ~ P",k)l(^^)  + 9i{p?+l,k  + P",k)l^  + 

*"[^(“r+l,*  + “"-I,*)  + S^^^lk+l  + “l!/b-l)l 
~^?ik  ~ (P?,k+1  ~ P?,k)ti^^)  + 9z{p"k+l  + P?,k)/^  + 

+ <-l,k)  + ^(“’"*+1  + <*-l)l 


We  define 


r^2uSt{^+±) 


Then,  for  the  second-order  time  steps  using  the  DuFort-Frankel  scheme 


= (1  + rr^i^StUPi,  + (1  - r)«rr'] 

= (1  + rr^[2StWp,  + (1  - r)«,.V‘] 

and  for  the  first-order  time  steps  (initially)  using  DuFort-Frankel 


= StWp,  + (1  - r)u;.’i 


4.3  Stability 

The  Linear  stability  of  each  of  the  two  schemes,  the  DuFort-Frankel  or 
Fromm  method  and  the  lagged-dissipation  method,  has  been  checked  for  a 
single  convection-diffusion  equation.  Derivation  of  the  dispersion  relation 
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in  each  case  is  straightforwcird;  however,  analysis  of  the  stabihty  (loca- 
tion of  the  roots  of  the  quadratic  in  the  dispersion  relation)  in  the  Fromm 
method  is  not  easy,  and  we  have  resorted  to  a theorem  derived  by  Miller 
[16]  referenced  in  the  book  of  Peyret  and  Taylor  [17]  to  guarantee  the 
stability.  Assessing  stabihty  and  study  of  the  dispersion  relation  for  the 
lagged-dissipation  method  is  much  simpler. 

We  start  with  a single  convection-diffusion  equation 


09  ^09  09 

+ 


eV^9 


where  9 is  the  single  dependent  variable,  U and  V are  the  loccd  hnearized 
vector  velocity  components,  e is  the  dissipation  coefficient  and  is  the 
Laplacian.  define 

d(t"yXj,Zk)  = 6% 

and  let  Xj  = j8x,  Zk  = kSz,  = nSt. 

First,  consider  the  DuFort-Frankel  or  Fromm  method,  and  difference 
the  convection-diffusion  equation  as  follows: 


9^V  - 


J,k 


J,k 


u 


w 


8z^ 


['•'(^"+1,*  + - (1  + '■"W  + = 0 


where  r = 8zl 8x.  Now  let 


9^3, k = A„exp[i(5z(;7(^Ax)  + 


Then 


“h  iAfi 


8z'^ 

2Q8t,Ur 
8z 


8z^ 

. ,8z  W . 8z. 


28i€,  o . 8z  ,8z, 

cos(— + cos(-] 


where  Q = \/V'^  -f  W'^, 
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Define 


28ie{l-^r^) 

6z^ 

^_2QU.Ur  . It  ,^W  . I, 


then 


An+i{l  + d)-  An-i{l  -d)-  An{Pd  - 2iT)  = 0 
or,  with  An  = and  dividing  through  by  9'^~^ 

S\1  4-  d)  - 0{Pd  - 2iT)  - (1  - d)  = 0 

Now,  according  to  the  theorem  of  Miller  [16],  if  we  define 

f(e)  = 9^-  9^^—^  - IzA 
' l + d l + d 

and  a polynomial,  constructed  from  /(^),  a<s  follows: 


m = 1 - ^ 


Pd  + 2iT  1 — d ^2 


l + d 


l + d 


0^ 


then,  we  can  construct  a polynomial  fi(0)  of  one  degree  lower  than  /(^), 
which  in  this  ca^e  is  of  first  degree,  as  follows: 

l.r. 


= 9 


1 

cs 

1 

rH 

r 

1 

1 

\^  + d)  \ 

l+d 

= ^imm  - /(o)/(^)i 

Pd(l  + - 2iT  fl 


l+d 


l + d 


The  theorem,  proved  in  Miller  [16],  states  that,  if  |/(0)|  = 1 > 1/(0)|  = 
(1  — d)/(l  + d)  and  |^i|  < 1,  where  9i  is  the  zero  of  /i(^i)  = 0,  then  both 
of  the  zeros  of  f{9)  = 0 are  such  that  |^|  < 1.  Now 


9i  = P/2  - iT 
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With  the  assumption  that  r = 1,  we  can  write  P and  T as  follows: 

2QSt  Sz  1 1 8z  11 

T = _ cos  a cos[  -( - - -)]  sm[  y ( ^ + j;)l 

_ ^ ,Sz . 1 1 ..  .Sz,  1 1 

p = 2cos[-(-+-)]cos[-(---)] 


If  the  Courant  condition,  < 1,  is  satisfied,  then  |^i|  < 1,  and  Fromm’s 
method  is  stable,  a conclusion  determined  by  MiUer  [16]  and  restated  by 
Peyret  and  Taylor  [17]  for  the  one-dimensional  convection-diffusion  equa- 
tion. 

For  the  lagged-diffusion  scheme,  a similar  analysis  can  be  performed; 
we  start  with  the  difference  form  of  the  convection-diffusion  equation: 


u 


w 


+ {ejil,  + erk-,)  - (i  + = o 

The  equation  for  the  amplitude  A„  in  this  case  becomes 


-'i-n+l[l  + 


28e{l  -f-  r^) 


8z‘^ 

2Q8t^Ur 
8z 


— ^n-l[l  — 


28e(l  H-  r^) 


8z^ 

^ . 8z  W 8z  . 

^ 28te,  2 / X 

[r  cos(— )+cos(-)] 


or, 

An+i(l  d)  — A„_i(l  — d)  — An-iPd  An2iT  = 0 

where  the  same  notion  as  for  the  Fromm  method  is  used.  Substitution  of 
An  = and  division  by  yields  the  dispersion  relation  for  the  lagged- 
diffusion  scheme 


0^1  + d)  + e2iT  - [1  - d(l  - P)]  = 0 

The  roots  of  this  quadratic  equation,  the  amplification  factors  for  the 
scheme,  are 

9 = Y^[-iT  ± v^-T2  + (1 + </)[! - </(!- P)]] 
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If  the  Courant  condition,  < 1,  is  satisfied  and  d = < 1,  then, 

for  both  roots, 


1 -d(l  -P) 
l + d 


< 


1 


since  |P|  < 1,  and  the  scheme  is  stable. 


4.4  The  Pressure  Equation 

The  incompressibility  condition  is 

{ui^k  - - Wi^k-\)ISz  = 0 (20) 

and  this  condition  can  be  used  to  derive  the  Hnear  algebraic  equation  system 
for  the  pressure.  Away  from  boundaries,  that  is  for  2 < i < I — 1^2  < k < 
— 1,  we  have 

Pi+l,k  ~ + Pi-1, k Pi,k+1  — '^Pi,k  + Pi,k-l  _ 

6'^X  8'^Z 

~‘r{~F'x,i,k  + Fx,i-i,k) / 8x  + {—Fz^i^k  H“  Fz,i^k-i)/8z 
+9xl—{pi+i,k  + Pi,k)  + (pi,k  + Pi-i,k)]/i2Sx)  + 

9z[  — {pi,k+l  + pi,k)  + {pi,k  + Pi,k-l)]/ {28  z)  + Tik 

Here  the  matrix  Tik  is  due  to  the  viscous  terms  (T  is  proportional  to  u), 
and,  as  above,  the  lack  of  a superscript  implies  that  all  quantities  are  to  be 
evaluated  at  time  level  n. 

At  the  boundaries,  we  have  four  edges  and  four  comers  at  which  to 
determine  the  equations.  Along  each  of  the  edges  and  corners,  we  system- 
atically write  out  the  equations  using  the  conditions  normal  to  the  boundary 
that  there  is  no  velocity  nor  velocity  flux  and  the  pressure  is  hydrostatic. 
For  1 < k < K 


uok  = 0 
ujk  = 0 
Fx,ok  = 0 
Fxjk  = 0 
Pok  = {8xgx/2){pik  + POAr) 
Pi+i,k  - Pik  = {8xgxl2){pi+i^k  4-  pik) 
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For  1 < i < I 


w,o  = 0 
WiK  = 0 
^z,iO  = 0 
Fz,tK  = 0 

Pil  - Pio  = {Szgz/2){p,i  -h  /5,o) 

Pi,K+l  — PiK  — + PiK) 

Also,  we  must  either  apply  free-slip  or  no-slip  BC.  For  0 < /?  < AT, 


and  for  0 < ^ < / 


h?/" 

“i,o  — ^“i,l 


u 


n 

i.K 


= ±U 


n 

i,K-l 


where  the  plus  sign  corresponds  to  free-slip  and  the  minus  sign  to  no-slip 
BC. 

To  account  for  tangential  boundary  conditions,  we  have  included  the 
matrix  T^.  The  expression  for  the  matrix  T,,j  is  straightforward  but  tedious 
to  derive.  In  the  paragraphs  below,  we  indicate  a derivation  for  this  matrix. 

Define  the  divergence  bls  follows 


D"k  = (Kk  - 


(21) 


If  we  impose  the  condition  on  all  interior  cells  (at  the  center  where  the 
divergence  is  defined)  that  the  divergence  be  zero  after  any  time  step,  then 
= 0 translates  into  the  condition  that 


{Ullk  - - W",_,)l6z  = 0 (22) 

where  and  are  the  right  hand  sides  of  the  momentum  equations 
defined  above.  Substitution  for  these  expressions  for  the  DuFort-Frankel 
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scheme  (the  expressions  for  the  lagged-dissipation  scheme  can  be  similarly 
derived)  then  yields 

Pi+i,k  — + Pi-1, k Pt,fc-t-i  - 4-  Pi,k-i  _ 

~^{~^x,i,k  + Fx^i^i^k)/Sx  + {—Fg^i^k  + Fz^i^k-l)/Sz 
~^9x[—{pi+i,k  4-  Pi,k)  + {pi,k  + Pi-\,k)]/i2Sx)  -h 
9z[  — {pi,k+l  H-  Pi,k)  + (pi,Jt  + P,-,A:-i)]/(2^z)  — 

■(“i.Jfc  "I"  ^i~2,k)  ~ “i-l,*-!)] 


(5^2 


+ “'■-i.-fc)  + 7^(“’a+i  + <k-i) 


^(“’"+1,*-!  + “'"-l,fc-l)  - -^(Kk  + <fc-2)l 


If  we  define 


D?.k  = {D?+i.k  + Di-i,k)li^'  + {D?m  + D",.,)/Sz^  (23) 

then,  the  pressure  equation  can  be  rewritten  els 

Pi+l,k  — *^Pi,k  H-  Pi~.l,k  Pi\fc+1  ~ ^Pi,k  + Pi,k-1  __ 

8'^X  8'^Z 

~^{~Fx,i,k  4-  Fx,i-\,k) / 8x  + {—Fz,i,k  4-  Fz^i^k-i)/ 8z 
-\-9x[~{pi+i,k  4-  Pi,k)  4-  {pi,k  4-  pi-i^k)]/i^8x)  -4 
9z[—{pi,k+i  4-  Pi,k)  4-  (pi,k  4-  pi,k-i)]/{^8z)  + i^D^^k 


Now,  in  cells  away  from  boundaries,  2 < i < I — 1 and  2 < A;  < A"  — 1, 
by  imposition,  aH  of  the  = 0 and  therefore  all  of  the  = 0.  Hence 

T-jt  = I'D'^k  — only  in  cells  adjacent  to  the  boundaries  (and  in 

corner  cells)  where  the  matrix  will  in  general  be  nonzero.  To  derive  an 
expression  for  this  matrix,  consider  the  cells  adjacent  to  the  left  boundary, 
i = 1,2  < k < K — 1.  Along  this  boimdary,  Uq^.  = 0 for  all  n so  that 
Uq  i^  = 0 implies 


1 


Subtracting  this  equation  from  the  general  equation  for  pressure  with  i = 1 
yields  the  equation 


P2,k^  Pi, ^ gl,fc+l 


Sx^ 

— 9x(P2,k  + Pi 


^ — \-^z,l,k  — JT  OZ 

^k)l{2Sx)  4-  gz[—(pi,k+i  4-  pi,k)  -f  {p\,k  4-  pi,k-i)]/i2Sz) 


Now,  we  impose  the  condition  that  the  divergence  be  zero  in  the  ficti- 
cious cell  outside  of  the  left  boundary  (which  we  are  free  to  do)  Dqj^  = 0. 
This  condition  implies  that 


= (<^'o.k  - “’o,fc-l)/^Z 


Similarly,  the  condition  that  the  divergence  be  zero  in  the  cell  adjacent  to 
the  left  boundary  = 0 implies  that 


u 


Now,  for  free-slip  BC, 


= -(Kk  - 

^0,k  — ^l,k 


and,  for  no-sHp  BC, 

^0,k  — ~ 

Using  all  of  these  relations  implies  that 


~ ^^l,k 

where  the  minus  sign  applies  for  free-sHp  BC  and  the  plus  sign  for  no-slip 
BC.  Using  the  imposed  condition  that  Dq^.  = 0 implies  that  = 0. 
Hence 

t'D?,k  + + “-!.<=)  = 

or,  = 0 for  free-sHp  BC,  while  j 8x^  for  no-slip  BC. 

In  a similar  fashion,  the  matrix  can  be  derived  at  all  other  cells 
adjacent  to  boundaries.  In  this  fashion,  we  find  that  for  free-slip  BC,  T,’J  = 0 
for  all  i and  k^l<i<I^l<k<K.  For  no-slip  BC,  the  matrix  is  still 
zero  at  cells  removed  from  the  boundary;  i.e.,  T,’J  = 0 for  2<z</  — 1,2< 
k < K — 1.  At  cells  adjacent  to  the  boundary,  i.e.,  edge  or  corner  cells,  the 
matrix  is  specified  below. 
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4.4.1  Edges 

For  A;  = 1, 2 < ^ < / — 1, 

Pi+1,1  — + Pt-1,1  Pi, 2 — Pi,i  __ 

Sx^  6z^  ^ 

+^x["-(pi+i,i  + Pi,i)  + (pi,i  + Pi=i,i)]/{25x) 

~'9z{pi,2  + Pt,i)/(2^z)  + r,-,i 

For  k ^ K,2  < i < J — 1,  we  have 

Pi+1,K  — ^Pi,K  + Pi~l,K  -^Pi,K  4-  Pi,K-l  _ 

Sx"^  Sz'^ 

4-(“F’x,i,ic  + Fx,i=i,K)l 8x  H-  Fz,i,K-\l 8z 

-^[-~9x{Pi-¥l,K  4-  pi,K)  + 9x{pi,K  4-  p,_i,^)]/(2^®) 
9z{Pi,K  + Pi,K~l)/{^8z)  -h  Ti^K 
Next  do  the  vertical  edges;  i = l,2<^<Ff--l. 

P2,k  - Pi,k  Pi,k+i  - 2pi,A;  4-pi,fe-i  __ 

8x'^  82^^ 

—Fx,i,k/ 8x  4-  {—Fz,\,k  4-  Fz,i,k--i)l 8z 

—9x{p2,k  4-  pi,k)/{28x)  4-  [— Pz(pi,jfc+i  4-  Pi,k)  4- 
9z{pi,k  4-  Pi,fc-i)]/ {282)  4-  Ti^k 

For  i = /,  2 < ib  < iF  - 1. 

—Ppk  4-  Pi-i,k  Pi,k+i  2p/,fc  4-  pi,k-i  __ 
8x^  82^ 

+Fx,i-i,k/8x  + {-’Fz,i,k  4-  Fz,i,k-\)l 8z 

+9x{pi,k  4*  pi-\,k)l{28x)  4-  [■~9z{pi,k+i  4-  pr,k)  4- 
9z{pi,k  4-  p/,*;_i)]/(2^2:)  4-  Ti^k 


4.4.2  Corners 

Finally,  do  the  comers.  First,  i = 1,  — 1 

P2,l  ~P\,i  , Pi, 2 “Pl,l  __ 

4- 


8x^  ' Sz^  ■F’-.u/'5* 

-Fz,i,il8z  - Px(P2,i  4-  p\,i)l{28x) 
—9z{p\a  4-  pi,i)l{28z)  4-  Ti,i 
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= Fxj-i,ilSx 


Next,  2 = /,  A;  = 1. 

~Pi,\  + P/-i,i  . PL2  — P/,1 
8x^  ^ 8z^ 

8z  + Px(p/,i  + pi~i,\) ! {28x) 

—9z(pi,2  H-  Pi,i)/{28z)  -f-  T/,1 

Next  i = l^k  = K . 

P2,K  - P\,K  , -Pl,K  + P\,K-l  „ 

6x2  + = -F.,..k/Sx 

-\-Fz^i^k-\/ 8z  — gx{p2,K  + Pi,k)/{28x) 

-^9z{pi,K  H"  P\,k-\)I{2‘8z)  -h  Ti^k 

For  i = I^k  = K. 

-PI,K  + Pl-1,K  , -Pl,K  + Pl,K-l  „ /X- 

+ s? = 

+^z,/,K-i/ 8z  + 9x{pi,K  + Pi-i,k)I  {28x) 
-\-9z{pi,K  + Pi,k-i)I{‘28z)  + T/,k 


Thus,  the  viscous  matrix  is  defined  as  follows:  For  free-slip  BC,  for  1 < i < 
1,1  <k<K 

Ta  = 0 

For  no-shp  BC,  for  2<2</  — l,2<A;<i^  — 1 


Ti,k  = 0 


for  2 < k < K — 1 

^ 2v 
“ 6x3“'’* 
_ 2u 


and  for  2 < i < / — 1 


.r,  2i/ 

T...  - 


Ti,K  = 


2u 


8z^ 
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Finally,  at  the  corners 


2i/ 

2u 

Ti.i 

= 

2u 

2u 

Ti^k  — 

1 

1 « 

+ 

Ti,i  = - 

2z/ 

2u 

Once  again,  no  superscript  implies  quantities  evaluated  at  time  level  n. 

5 Results 

In  this  section  samples  of  results  generated  using  a code  implementing  the 
algorithm  described  above  are  presented  and  discussed.  The  results  repre- 
sent some  problems  of  interest  to  scientists  concerned  with  enclosure  fires. 
The  aim  of  the  computations  is  to  investigate  basic  phenomena  which  occur 
in  enclosure  fires  using  current  computational  and  graphical  capabilities. 
The  computations  were  performed  on  the  Convex  C120  in  the  Computa- 
tional Combustion  Facility,  which  is  a joint  facility  of  the  Computing  and 
Applied  Mathematics  Laboratory  and  Building  and  Fire  Research  Lab- 
oratory at  NIST.  The  graphics  were  generated  on  four  Silicon  Graphics 
Personal  Iris  4D20s,  which  are  also  part  of  this  facility. 

The  scale  of  these  computations  is  substantial.  Runs  have  now  been 
performed  using  almost  one-half  million  grid  points;  they  require  nearly 
50  megabytes  and  take  up  to  24  hours  of  CPU  time.  They  are  performed 
in  64  bit  arithmetic.  If  we  were  to  save  any  major  fraction  of  the  data 
generated,  we  would  be  overwhelmed.  Therefore,  a significant  effort  has 
been  expended  in  trying  to  select  only  those  data  required  to  understand 
the  phenomena  being  studied.  We  have  concluded  that  Lagrangian  particle 
tracking  for  transient  phenomena  is  the  most  convenient  method  for  saving 
and  examining  data  from  our  computations,  and  this  procedure  is  described 
below. 

The  flows  are  visualized  by  introducing  Lagrangian  particles  into  the 
heat  source  representing  the  fire.  These  particles  are  introduced  by  means 
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of  a random  number  generator  which  mimics  the  spatial  distribution  of 
heat  in  the  source.  Therefore,  in  general  in  the  results  shown,  the  particles 
are  introduced  according  to  a Gaussian  distribution  centered  at  the  center 
of  the  heat  source  and  with  a half-width  equal  to  that  of  the  source.  The 
particles  are  introduced  at  regular  intervals  of  time,  again  mimicing  the 
temporal  dependence  of  the  heat  source,  which  is  “turned  on”  rather  quickly 
(on  the  dimensionless  time  of  the  computations)  to  a constant  rate  of  heat 
release.  Each  set  of  particles  is  then  allowed  to  move  with  the  flow  field 
so  that  one  can  see  the  movement  of  the  “smoke  and  hot  gases”.  The 
graphics  are  generated  in  color  with  the  color  representing  the  temperature 
of  the  particle.  In  this  fashion,  one  can  also  discern  where  the  flow  is 
hot  and  where  it  has  cooled  by  mixing  and  conduction.  The  color  and  the 
dynamic  aspect  of  the  graphics  (a  sequence  of  frames  such  as  those  exhibited 
below  is  displayed  rapidly  on  the  computer  screen)  enhance  enormously  the 
information  which  can  be  extracted  from  the  computations. 

The  amount  of  computer  memory,  the  CPU  time  and  the  data  storage 
all  increase  as  the  resolution  of  the  computations  increases.  If  /,  K are  the 
number  of  grid  points  in  the  x and  y directions,  then  the  memory  require- 
ments scale  B.S  IK.  If  we  assume  that  I = K,  then  the  CPU  time  scales 
as  and  the  disk  space  scales  as  also  if  we  retain  data  computed 
on  the  grid.  If,  on  the  other  hand,  we  retain  particle  data,  then  the  disk 
requirements  can  be  decoupled  from  the  grid  size  and  are  dependent  on  the 
number  of  Lagrangian  particles  introduced  and  the  number  of  times  these 
particle  data  are  dumped.  (Some  additional  time  is  required  to  compute 
the  particle  trajectories;  however,  as  long  as  the  number  of  particles  intro- 
duced is  small  compared  to  the  number  of  grid  points,  this  additional  CPU 
time  is  relatively  small.  In  fact,  for  most  computations  we  have  done,  this 
additional  time  amoimts  to  only  ten  percent  or  less  of  the  total  CPU  time.) 

The  resolution  of  the  computation  determines  the  Reynolds  number  of 
the  flow  which  can  be  calculated,  and  this  is  a very  important  feature  of 
the  computations.  Since  the  size  of  the  Reynolds  number  for  the  2-D  flow 
scales  as  with  the  grids  that  we  have  been  able  to  use,  approximately 
one-half  million  grid  cells,  flows  with  a Reynolds  number  of  nearly  one- half 
million  can  in  principle  be  computed.  We  have  not  yet  computed  flows  with 
this  large  a value  of  the  Reynolds  number,  primarily  because  we  have  been 
examining  non-square  geometries.  However,  we  have  computed  results  for 
flows  with  Reynolds  numbers  of  up  to  1 x 10^.  One  comforting  aspect  of  the 
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algorithm  and  the  code  is  that  it  will  “bomb”  if  the  Reynolds  number  is  too 
large  for  the  resolution  of  the  grid.  This  feature,  plus  the  fact  that  there 
are  no  adjustable  parameters  in  the  code,  lead  us  to  have  great  confidence 
in  the  predictive  capability  of  the  computations. 

5.1  Stratified  Background  Fluid 

A problem  of  some  general  interest  is  that  of  fire  in  an  enclosure  where  the 
background  fluid  has  been  stratified  [9],  [10],  [11],  [12].  Stratification  occurs 
naturally  in  any  room  that  is  relatively  quiescent.  In  the  summer,  solar 
heating  generally  produces  significant  stratification,  and  during  the  winter, 
unless  the  stratification  is  disrupted,  normal  house  heating  also  produces  it. 
A small  fire  under  this  circumstance  will  have  its  induced  flows  significantly 
altered  by  the  ambient  stratification.  Also,  if  a fire  has  started,  produced  a 
stratified  upper  layer  in  a room,  temporarily  died  down,  and  then  revives 
again,  the  plmne  dynamics  wiU  be  affected  by  the  stratification. 

The  first  example  is  of  a heat  source  (a  line  source  in  these  2-D  calcula- 
tions) at  the  floor  of  a rectangular  enclosure.  The  aspect  ratio  of  the  room 
is  0.8,  the  stratified  layer  occurs  at  0.6  the  height  of  the  enclosure  and  the 
difference  in  density  between  the  lower  and  upper  layers  is  5 times  the  char- 
acteristic density  induced  by  the  heat  source.  The  heat  source  is  centered 
on  the  floor,  one-half  unit  (the  length  unit  being  the  height  of  the  room) 
from  the  left  wall  with  a Gaussian  half-width  of  approrimately  one-tenth 
unit.  Four  frames  of  a sequence  are  shown  in  Figure  1.  The  first  frame 
shows  the  starting  plume  at  2.6  dimensionless  time  units.  Note  both  the 
symmetry  of  the  starting  plume  at  this  time  and  the  detailed  structure  of 
the  head  of  the  starting  plume.  This  structure  is  due  to  the  high  resolution 
of  the  computation,  which  has  600  x 480  grid  cells,  and  has  a Reynolds  num- 
ber of  1 X 10®.  The  second  frame  in  Figure  1 shows  the  plume  at  3.2  time 
units.  At  this  time,  the  plume  has  encountered  the  upper  stratified  layer; 
part  of  the  plume  penetrates  the  layer,  but  most  of  the  plume  is  deflected 
laterally.  Note  the  asymmetry  which  has  begun  to  appear  due  to  the  fact 
that  the  heat  source  is  not  symmetrically  placed  in  the  room.  In  the  third 
and  fourth  frames  the  plume  continues  to  be  deflected  laterally  because  the 
upper  layer  is  so  strongly  stratified.  Again,  note  the  asymmetry  in  the  flow 
field. 

Computations  have  also  been  carried  out  for  other  levels  of  stratification 
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of  the  upper  layer.  As  this  stratification  becomes  stronger,  the  amount  of 
penetration  of  the  plume  into  the  upper  layer  is  reduced.  On  the  other 
hand,  as  this  stratification  becomes  smaller,  more  penetration  occurs  with 
the  plume  reaching  the  ceiling,  and  finally,  with  the  stratification  having 
little  effect  when  it  becomes  very  weak. 

5.2  Gravity  Currents 

When  lighter  fluid  is  introduced  at  the  top  of  an  ambient  fluid  (or  heavier 
fluid  at  the  bottom  of  ambient  fluid),  the  light  fluid  spreads  over  the  am- 
bient. The  process  by  which  this  spreading  takes  place  is  called  a gravity 
current.  Gravity  currents  are  well-known  in  the  geophysical  literature  and 
have  been  studied  both  experimentally  and  theoretically  for  many  years 
[13].  We  have  used  our  2-D  model,  since  gravity  currents  are  primarily 
two-dimensional  phenomena,  to  compute  gravity  currents  with  very  high 
resolution  and  relatively  high  Reynolds  numbers  (up  to  3 x 10"^).  In  Fig- 
ure 2,  we  show  several  frames  from  a sequence  of  computed  frames  of  a 
gravity  current  generated  at  the  left  of  an  enclosure  and  progressing  to 
the  right.  The  gravity  current  has  the  weU-known  characteristic  “head” 
and  progresses  at  a nearly  uniform  speed.  The  enclosure  is  11  units  long 
and  one  unit  high,  and  has  1584  cells  in  the  horizontal  direction  and  144 
in  the  vertical  (228,000  grid  cells).  The  Reynolds  for  this  computation  is 
1 X lO'^.  From  these  computations,  we  can  compare  results  with  experimen- 
tal measurements  carried  out  elsewhere.  First  we  can  compare  the  velocity 
of  the  head  of  the  gravity  current  and  determine  its  variation  with  Reynolds 
number  for  example  and  later  we  may  be  able  to  make  more  detailed  com- 
parisons, comparing  velocity  or  density  profiles  through  the  gravity  current 
for  example. 

5.3  The  ^‘Trench  Effect” 

Fires  in  buildings  involve  the  transport  of  heat  and  mass  by  gravity-induced 
or  buoyant  convection.  Generally,  this  convection  occurs  in  rectangular 
enclosures  where  the  direction  of  gravity  is  parallel  to  the  surfaces  of  the 
enclosure,  the  wails.  However,  under  certain  circumstances,  such  as  a fire  in 
a stair  well  or  an  escalator,  the  enclosure  may  be  sloped  relative  to  gravity. 
A very  important  example  of  a fire  in  a sloped  enclosure  was  the  devastat- 
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ing  fire  in  the  King’s  Cross  underground  station  in  England  in  1987,  where 
there  was  significant  loss  of  fife  as  well  as  property  damage.  Numerical 
simulation  of  this  fire  uncovered  an  unexpected  phenomenon  which  caused 
a very  rapid  spread  of  the  fire  and  led  to  much  of  the  devastation  [14], 
[15].  This  phenomenon  was  termed  the  “trench  effect”,  and  caused  some 
controvery  during  investigations  of  the  King’s  Cross  fire  in  England.  The 
phenomenon  was  ultimately  confirmed  by  experiments  and  additional  sim- 
ulation, but  transient  aspects  of  the  fire  simulation  are  stiU  of  interest.  We 
have  examined  this  phenomenon  using  our  2-D  code  for  Reynolds  numbers 
of  interest  for  enclosure  fires. 

In  Figure  3,  we  present  a sequence  of  frames  from  a 2-D,  high-resolution 
computation  for  the  flow  generated  in  an  enclosure  by  a heat  source  (fire) 
located  on  the  floor  one  quarter  of  the  length  form  the  left  wall.  The  resolu- 
tion is  1296  X 324  (419,904  grid  cells)  in  this  4x1  enclosure.  The  Reynolds 
number  for  this  computation  is  5 x 10"*,  and  free-slip  and  adiabatic  boimd- 
ary  conditions  have  been  imposed.  The  enclosure  is  inclined  35  degrees 
with  respect  to  horizontal.  In  this  flgure,  the  plume  rises,  but  is  bent  back 
toward  the  back  wadi.  After  the  hot  gases  hit  the  ceiling,  they  progress  both 
toward  the  back  wall  and  up  the  ceiling  toward  the  high  end.  However,  the 
hot  gases  leaving  the  heat  source  are  pinned  along  the  floor  and  form  a 
hot  gas  jet  which  progresses  up  along  the  floor,  shedding  hot  gases  near  its 
front;  this  phenomenon  we  interpret  as  the  “trench  effect”.  These  results 
were  imexpected  and  were  the  reason  that  the  “trench  effect”  caused  so 
much  controversary. 
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A Appendix:  Particle  Tracking  Algorithm 

A.l  Particle  Injection 

The  initial  location  of  each  particle  is  selected  at  random  from  a distri- 
bution function  which  is  normal  in  the  x-direction  with  mean  value  equal 
to  the  location  of  the  center  of  the  heat  source  and  variance  equal  to  the 
mean  width  of  the  heat  source.  In  the  z-direction  the  distribution  func- 
tion is  either  normal  with  characteristics  similar  to  that  in  the  x-direction 
or  exponential.  In  all  cases,  the  distribution  functions  mimic  the  spatial 
distribution  of  the  heat  source.  Also,  particles  are  injected  at  a rate  that 
follows  the  rate  of  heat  addition. 

A. 2 Interpolation  of  the  Fields  at  Particles  Locations 

Take  the  2-D  case  and  let  a?”  and  Zj  be  the  coodinates  of  particle  < 
j < Np  at  time  t”.  Let  the  grid  be  our  standard  uniform  grid, 

Xi  = iSx,  0 < i < I 

Zk  = kSz^  0 < k < K (24) 

Let 

I,  = [x^/Sx  + 1/2] 

ICj  = [z’^/Sz  + 1/2] 

ri  = x1ISx-{I,-\l2)  (25) 

3^  = z"l8z  - {ICj  - 1/2) 

and  let 

r.^lx^/sx] 

K.r  = [zysz] 

r'  = x"/5z  - I,  (26) 

s'j  = Zj  /Sz  — Kj 

where  [...]  is  the  integer  part  of  the  quantity  in  square  brackets. 
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Define  the  function  F as  follows: 

^(-*^1,1 5 -*^1,2 1 -*12,11  ^2,2,  ^)  = ^4i,i(l  — r)(l— ,s)H- .42,1  r(l—s)+Ai, 2(1  — r)^+A2, 2^5 

and  let  Uj,  iDj,  pj  be  the  velocity  components  and  the  density  respectively 
at  the  location  of  the  jth  particle.  Then  by  biHnear  interpolation: 


w 


~n  IT /'ll”  11  ” 11”  11”  -n'  a\ 

Uj  — -r  Uj;^2,a:j+ii  ^ij+i,x:j+2i  “jj+2,;Cj+2i 

n ZT'/'iii”  im”  11,”  11,”  T,  a'\ 

j — r \^2j  + i,lC'j-\-l^^Ij+2,K'j  + l'^Ij  + l,/C'j+2^  ^Jj+2,K'j+2^  ^ ) 

~n  77/  .n  n n n / f\ 

Pj  — ^ lPjj+i,x:'.+ii  Pi'-i-2,x:'.+i’  Pjj+i,/c'+2i  Pi'+2,x:'.+2i  ^ 1 ; 

A. 3 Time  Advancement  Scheme  for  Particles 

Now,  each  particle  is  advanced  through  the  differential  equation: 

dxj 


(27) 

(28) 
(29) 


dt 

dzj 

dt 


~ ^3) 


(30) 


The  particles  are  advanced  in  time  according  to  a second  order  Runge- 
Kutta  time  integration  scheme  as  follows. 


Then,  with 


rj  = 
3 1 — 


hl^  = u"(x",  z")St 
K,i  = w'j{x’j,Zj)St 


I,  = [(x’;  + hl^)/6x  + l/2] 
K,  = [(2"  + hl,)l6z  + 1/2] 
{x"  + hl,)ISx  - {I,  - 1/2) 
{z]  + hl^ysz  - (IC,  - 1/2) 


(31) 


and  let 


i;.  = [(x"  + hiy/sx] 
)C’j  = [(2"  + Ky/5z] 
r'  = (x"  + h",)/Sx  - Ij 
= iz"  + hiy/Sz  - Ki 


(32) 
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The  new  extrapolated  velocity  components  at  location  H-  h"  + h” 

3 ^i3'  J z,j 

are 


~ “i;+2,AC>+l’ “l;+l,Kj+2)  “rj+2,/Cj+2i 

~ “’rj+2,/cj+i!  “^ij+i,;c^+25  “’Ij+2,a:'+2> 


r',^) 

r,  s') 


Finally,  the  location  of  the  particle  at  time  level  n + 1 is 

x"+'  = x]  + 0.5(uj"5<  + /!"_,) 

= z"  + 0.5(w’]St  + hlj) 


(33) 

(34) 

(35) 


(36) 


The  accuracy  of  this  particle-tracking  routine  was  tested  using  a steady- 
state  flow  with  vorticity  in  a rectangular  enclosure.  This  flow  field  was 
originally  derived  to  test  the  algorithms  that  integrate  the  fluid  equations 
and  is  described  in  [18].  In  the  flow,  particles  follow  streamlines  which 
close  upon  themselves,  a sensitive  test  of  the  quality  of  the  particle-tracking 
integration  scheme. 
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Figur6  1 Buoyant  Convection  Induced  by  a Heat  Source  in  an 

Enclosure  with  Background  Stratified  Fluid 


Frame  2-D,  frame  no.  5,  time  = 5.0000 


Frame  2-D,  frame  no.  10,  time  = 10.0000 


Frame  2-D,  frame  no.  20,  time  = 20.0000 


Frame  2-D,  frame  no.  30,  time  = 30.0000 


Frame  2-D,  frame  no.  40,  time  = 40.0000 


Figure  2 

A Gravity  Current 
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Pi  _ Buoyant  Flow  Induced  by  a Heat  Source  in  an 

yur0  o Enclosure  Inclined  at  35  Degrees 
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